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Chemical reaction rates must increasingly be determined in systems that evolve under the control of external 
stimuli. In these systems, when a reactant population is induced to cross an energy barrier through forcing 
from a temporally varying external field, the transition state that the reaction must pass through during 
the transformation from reactant to product is no longer a fixed geometric structure, but is instead time- 
dependent. For a periodically forced model reaction, we develop a recrossing-free dividing surface that is 
attached to a transition state trajectory [T. Bartsch, R. Hernandez, and T. Uzer, Phys. Rev. Lett. 95, 
058301 (2005)]. We have previously shown that for single-mode sinusoidal driving, the stability of the time- 
varying transition state directly determines the reaction rate [G. T. Craven, T. Bartsch, and R. Hernandez, 
J. Chem. Phys. 141, 041106 (2014)]. Here, we extend our previous work to the case of multi-mode driving 
waveforms. Excellent agreement is observed between the rates predicted by stability analysis and rates 
obtained through numerical calculation of the reactive flux. We also show that the optimal dividing surface 
and the resulting reaction rate for a reactive system driven by weak thermal noise can be approximated well 
using the transition state geometry of the underlying deterministic system. This agreement persists as long 
as the thermal driving strength is less than the order of that of the periodic driving. The power of this result 
is its simplicity. The surprising accuracy of the time-dependent noise-free geometry for obtaining transition 
state theory rates in chemical reactions driven by periodic helds reveals the dynamics without requiring the 
cost of brute-force calculations. 


I. INTRODUCTION 

Optimal control of reaction pathways in systems un¬ 
dergoing configurational changes can be achieved through 
forcing from tailored external fields. These fields can 
be tuned to induce specihc deformations on a poten¬ 
tial energy surface, providing control of state-to-state 
transitionsii^— In these processes, a formally exact classi¬ 
cal rate calculation can be obtained through modern-day 
transition state theory (TST)4i^— The principal assump¬ 
tions of TST are that (1) the distribution of energy states 
in the reactant configuration, and at the TS, are given 
by equilibrium distributions and (2) there exists a hyper¬ 
surface between reactant and product confirmations that 
is crossed only once by reactive trajectories during the 
traversal of a free energy barrier separating these basins. 
The TST reaction rate is calculated from the flux through 
this dividing surface (DS). If the DS is recrossed by re¬ 
active trajectories, TST will give an overestimate to the 
classical reaction rate. Determination of a recrossing-free 
DS therefore leads to classically exact TST rates systems 
in which equilibrium statistical mechanics is applicable. 

A phase space DS that is free of recrossings can be con¬ 
structed in conservative systems at energies close to the 
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reaction threshold. In systems with two degrees of free¬ 
dom, the optimal DS is the configuration space projec¬ 
tion of an unstable periodic orbit (PO)ii2r— In systems 
with higher dimensionality, the generalization of this PO 
is a normally hyperbolic invariant manifold (NHIM)fii“— 
The NHIM bounds the TS, being one less in dimension/^ 
It defines a recrossing-free surface at energies below bifur¬ 
cation thresholdsReactive trajectories are mediated 
by stable and unstable manifolds (reaction pathways) at¬ 
tached to the NHIM. These pathways persist even in re¬ 
actions whose state-to-state transitions are not dictated 
by purely configurational changes<^ 

In systems subjected to time-varying external forc¬ 
ing, the characterization of the NHIM as a hypersphere 
of constant energy breaks down. For example, field- 
matter interactions constitute processes in which energy 
is exchanged with a reacting system. These interactions 
lead to emergent and controllable behavior in assembly 
phenomena)^^— organic synthesis?^ protein folding;^ 
the detection of DNAf^ and photodissociation<^ Knowl¬ 
edge of the mechanism by which these interactions medi¬ 
ate reactive flow provides a methodological tool in the de¬ 
sign of molecular devices with unique functionalityi^^— 

Materials that undergo conformational changes in re¬ 
sponse to an external trigger offer examples of such emer¬ 
gent technology i^^i^^~ — Stimuli such as thermal varia¬ 
tions, electric helds, and photoinduction have been used 
as triggers for the conversion of chemical energy into me- 
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chanical worki^2i^“— Assemblies that convert chemical 
energy into directional motion can achieved through iso¬ 
merization reactions which are induced either from light 
or applied electric fieldsi ^^d^d^ In these responsive mate¬ 
rials, controlling the rate and pathway at which reactants 
transform to products is fundamental to harness mechan¬ 
ical actions for applicative purposes. 

The aim of this paper is to develop a rate theory for re¬ 
actions that are driven by periodic external fields in weak 
thermal environments. In the absence of noise, a dissipa¬ 
tive system that is periodically driven admits a DS that 
is free of recrossings4i This structure differs from the 
canonical view of the TS wherein the TS is a structure 
fixed in time at a saddle point on the potential energy 
surface. Here, we develop a rate theory based on reac¬ 
tive flux through this recrossing-free DS and the stability 
of the corresponding TS. In Ref. |4^ we found that the 
stability of the moving TS directly determines the reac¬ 
tion rate for single-mode sinusoidal driving. In conserva¬ 
tive systems, stability analysis is known to characterize 
molecular motions as well as determine the rate of con¬ 
figurational transitionsBuilding on our previous 
work, we test the viability of stability analysis to de¬ 
termine reaction rates in systems driven by multi-mode 
waveforms with no thermal driving. The extent to which 
the accuracy of the rate theory relying on the noise-free 
geometry persists in systems that are coupled to a ther¬ 
mal bath is also verified through inclusion and variation 
of the thermal driving strength. 

This outline of this paper is as follows: In Sec. m a dy¬ 
namical system is introduced to model barrier crossings 
in chemical reactions forced by periodic external fields. 
In Sec. imi a dividing surface that is recrossing-free is con¬ 
structed for this model in the absence of thermal driv¬ 
ing. Section IIVI contains analytical theories to predict 
the reaction rates of driven reactions by calculation of 
the reactive flux through this dividing surface for both 
globally non-linear and locally linear dynamics. Com¬ 
parison to the computational rates, computed from nu¬ 
merical integration of large ensembles of trajectories, is 
presented in Sec. |Vl Although not considered earlier, the 
effect of noise on the rate of these driven systems have 
also been addressed in Sec. ED We find that the rates 
computed from the noise-free geometry are accurate up 
to relatively large values of the friction and sometimes 
even in the thermal regime. 


II. MODEL DETAILS 

The interaction of an external field with a reactant 
species can strongly influence the mechanism and rate of 
a reactioni^i^i^ As a paradigmatic example of a chem¬ 
ical reaction driven under kinetic control, we consider a 
particle of unit mass moving along a reaction coordinate 
X. The trajectory of the particle begins at a position xq 
on the reactant side of an energy barrier that is moving 
in space under the influence of a time-dependent exter¬ 


nal field B(t). The chosen potential surface is the quartic 
form 

U(x) = -^u;i(x - E{t)f - \e{x - E{t))\ (1) 

The time dependent, instantaneous position of the mov¬ 
ing barrier top (BT) is specified by E{t). 

With the inclusion of additional non-conservative dis¬ 
sipation as well as stochastic driving forces, a particle at 
a phase space point F = (x, v) moving according to the 
potential o can be described by the Langevin equations 
of motion 

X = V, 

V = —7V -I- uj^(x — E{t)) + e{x — E{t))^ + 

( 2 ) 

where 7 > 0 is a dissipation parameter, Wb is the bar¬ 
rier frequency, and e is an anharmonic coefficient. Thus, 
for e 7 ^ 0 the coordinate of the particle is non-linearly 
coupled to the moving barrier. By restricting the anhar¬ 
monic coefficient to values e > 0 , there is a single max¬ 
imum in the potential located at the BT. The random 
fluctuating force ^a{t) is Gaussian white noise obeying 
the statistical properties 

{^a{t))=0 and (3) 

where a denotes a specific noise sequence. The strength 
of the noise is varied through the parameter cr > 0 . 

Depending on the geometry of C/(x), initial conditions, 
as well as the specific realization of the thermal environ¬ 
ment and the external field, a trajectory will either sur¬ 
mount the energy barrier and form product or remain on 
the reactant side. By calculation of the normalized flux 
of reactive trajectories through the the TS, the classical 
reaction rate for a system evolving through ([5]) can be 
obtained 1 ^ 

In this article, we consider periodic external driving of 
the form 

E{t) = al[ sin(a;t-I-(/?). (4) 

where fig C M is a finite set of frequencies. The wave¬ 
forms consist of a fundamental frequency D, and con¬ 
volutions of this fundamental with higher order partial 
frequencies. Three frequency sets fts are considered: the 
single fundamental frequency Jli = {D}, the fundamen¬ 
tal and the second partial frequencies fl 2 = {D, 2D}, and 
the fundamental, second, and third partial frequencies 
Da = {D, 2D, 3D}. The fundamental driving frequency 
D/ is D for Di and ^ 2 , and 2D for Da. The products in 
Eq. (|4|) for the three sets can be written as finite sums 

Ei{t) = asin(Dt -|- (f), 

E 2 {t) = -(cos(Dt)-I-cos(3Df-I- 2 (/ 3 )), 
a 

Eait) = ^(sin(2Dt -|- (p) -|- sin(4Dt + ip) 

+ sin( 6 Dt + 3ip) + sin(<p)). 
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FIG. 1. Functional forms of periodic driving E{t) for the single f2i (left), fl 2 (middle), and (right) are shown in the top 
row. The corresponding TS trajectory for each frequency set and anharmonic parameter values e G {0, 4, 8} are shown in the 
bottom row. Various extrema of E{t) are denoted by circles with arguments shown as dashed vertical lines. The corresponding 
extrema of x^{t) are denoted by circles and colored according to the respective e value. Units in all panels dimensionless. 


where the leading order terms exhibit the characteristic 
fundamental frequency. The maximum amplitude of each 
waveform is set to unity by adjusting the value of the 
parameter a accordingly. For the Hi, ^ 2 , and Ha sets, 
a = 1, a Ri 1.299, and a ~ 1.822, respectively. The 
functional forms of (O are shown in Fig. [T] 


III. THE TRANSITION STATE TRAJECTORY 


it is a moving saddle point to which stable and unstable 
manifolds can be attached. All trajectories that exponen¬ 
tially approach the TS trajectory as t —>■ oo are contained 
on the stable manifold. These trajectories will never de¬ 
scend from the BT region and therefore separate reactive 
from nonreactive trajectories in phase space. The unsta¬ 
ble manifold is formed from trajectories that approach 
the TS trajectory as t ^ —oo. The role of the unstable 
manifold is less important for the purposes considered 
here. 


The construction and existence of a structure whose 
configuration space projection is free of recrossings is de¬ 
pendent on the mechanism and geometry of a given reac¬ 
tion. For example, Mullen et have proposed that 

for ion-pair dissociation a no-recrossings DS does not ex¬ 
ist. This is in contradiction to earlier work by Truhlar 
and Garrett^ who proposed that through variation of a 
DS into an optimized orientation, recrossings could be 
eliminated. We have previously shown that for a period¬ 
ically driven system with no thermal driving, an optimal 
(recrossing-free) DS can be readily obtained. It is asso¬ 
ciated with an unstable PO in the region of the BT4i 
Moreover, a DS that is free of recrossings is known to ex¬ 
ist in thermally driven systems for the case of a harmonic 
barrier The time evolution of the configuration space 
projection of this DS has been termed the transition state 
trajectory4ii^i^“— It has not yet been proven is this is 
the only object which is free of recrossings over an arbi¬ 
trary finite time interval. Nevertheless, all that is needed 
here is its existence and the configuration projection of 
the TS trajectory defines a DS that is recrossing-free. 

The TS trajectory is a specihc trajectory that never 
descends into either the product or reactant regions, re¬ 
maining bounded to BT for all time. For the system ([2) , 


For an arbitrary driving E{t) of a harmonic (e = 0) 
potential, the equations of motion can be solved exactly 
and an exact form of the TS trajectory can be obtained. 
The eigenvalues of m 


Ku = , (6) 


correspond to the stable and unstable manifolds. The S 
functionals^i^ 




/ OO 

( 7 (r) exp(^(t — r)) dr : Re/r>0, 
+ / g^r) — t)) dr : Re/i<0, 

< J —OO 

( 7 ) 


obtained as a Green’s function solution, suppress the 
transient exponential factor in the solution and return 
only the equilibrium portion. In the absence of thermal 
driving (tr = 0), the TS trajectory for a harmonic barrier 
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FIG. 2. Phase space plots for a swarm of trajectories following the equations of motion © with various driving frequency sets 
and parameter values: f2i with 17/ = 5 and e = 1 for (a) cr = 0 (noiseless) and (b) a = 0.025 (thermal) driving, fls with f2/ = 4 
and e = 3 for (c) cr = 0 and (d) cr = 0.025. The initial position for every trajectory, xo = —0.1, is shown as a vertical solid 
white line. Reactive trajectories are colored in cyan and nonreactive trajectories are colored in orange. The TS trajectory T* 
is shown as a solid white curve. The critical velocity is indicated by a circle at the intersection of the dashed horizontal line 
and the line of initial conditions. The solid red curve is the critical curve and the solid blue curve is the harmonic critical 
curve Vc\e=o- Parameters for all panels are 7=1 and :/> = 0 in dimensionless units. 


can therefore be expressed 

xHt) = S[Xn.E-t\), 

vHt) = (As^[As, E-1] - A,5[Au, E- 1]). 


For the case of T-periodic motion of a harmonic bar¬ 
rier, the TS trajectory can be identified more easily by 
looking for a bounded solution to the equations of mo¬ 
tion. For the single frequency Hi case, the ansatz 

a;J(t) = Asin(r7t + (/?)-I- S cos(Ht -I-(/?), (9) 

yields the solution 


A = Ai, B = Bi, (10) 


where 


Ak = 


u;^(u;i + (km 
+ {uj^ + {kV,y 


Bk = a- 




At = a 


i + nytul’ 


Bi = 0 . 


( 11 ) 


(qfcH)^ + (w^ -b (k^YY 
In the absence of friction, 7 = 0, this simplifies to 

1 


( 12 ) 


In this case, the TS trajectory will oscillate in phase with 
the barrier, but with smaller amplitude Ai < a. 

For the r 22 and H 3 cases, ansatze can be constructed 
through Fourier series expansion of Eq. (|4]) yielding the 
solutions 


xlY) 


1(^1 cos(Ht) — Bi sin(Ht) 

— A 3 cos{3flt + 2 (p) + B 3 sin(3r77 + 2 ip)), 


( 13 ) 


and 

xl{t) = 7(^2 sin( 2 Ht + (fi) + B 2 cos(2r7t + (p) 

+ A 4 sin(4r2t + ip) + B 4 cos(4f2t + p) 

(14) 

— Aq sin(6f7t + 373 ) — Bq cos(60t -b 3:^3) 

-b asin((b>)), 
respectively. 

For periodically driven anharmonic barriers (e Y 0)j 
the TS trajectory rermains an unstable PO^ but it does 
not admit to an exact solution of the system of equa¬ 
tions Nevertheless, the phase space vector of the 
TS trajectory pl = (a;l(t),ul(t)) is a bounded solution 
to the equations of motion. To find this bounded solu¬ 
tion to arbitrary accuracy, numerical Newton-Raphson 
root finding methods were applied. 

The dynamics of x^{t), shown in Fig. [1] illustrate the 
result that the instantaneous position of the TS trajec¬ 
tory does not correspond to the energetic maximum of 
the potential surface. For dissipative systems (7 Y 0)j 
x^{t) will either lag behind in phase, as is the case for 
both the Hi and H 3 sets, or advance in phase as is the 
case for the H 2 set, with respect to motion defined by 
E{t). Also note that x^[t) oscillates with a smaller am¬ 
plitude than E{t). Thus, even for in-phase oscillations, 
e.g., when 7 = 0 , it will not correspond to the location of 
an energetic saddle point. Figure[T]also shows the depen¬ 
dence of x^{t) on the anharmonic parameter e. As e is 
increased, the curvature of the energy barrier increases. 
Non-intuitively, this results in a larger amplitude of oscil¬ 
lation for x^{t) to remain bounded to the BT. This trend 
persists for all H^. 

For dynamical analysis, it is advantageous to introduce 
a coordinate system which has a fixed point at the origin. 
In relative coordinates 

/S.X = x — x^ it), Av = v — v^{t), 


(15) 
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the equations of motion read 
Ai: = Av, 

+ + ( 16 ) 

Av = —jAv — U'{Ax + x^{t)) + U'{x^{t)). 

The relative equations of motion have a fixed point AT* 
at Ax = Av = 0, i.e., on the TS trajectory, and the sur¬ 
rounding vector field itself will now oscillate with period 
T, the same period as the driving. The TS trajectory has 
both a stable and an unstable manifold attached. In rel¬ 
ative coordinates, the directions of these manifolds will 
depend on time. 


IV. REACTION RATE THEORY 

In the TST formalism, the rate of a chemical reaction 
is given by the time-dependence of the conversion process 
from reactant to product (R —>■ P) where a DS in either 
configuration space or phase space separates the reactive 
constituents. The reaction rate can be obtained from the 
dynamics of the normalized reactive population (Pr —>■ 
Pp) either through analytical propagation of the phase 
space density of initial conditions or by treating large 
numbers of trajectories as discrete sets, and integrating 
the equations of motion. 

Consider a set of trajectories evolving through (l2|) that 
all have initial positions xq < a:^( 0 ) on the reactant side 
of the moving surface. The initial position distribution 
at time t = 0 is 5{x — xq) and the initial phase space 
density is 

Po{x,v) = 5{x - xo)q{v) (17) 

where q(y) is a Boltzmann distribution. The initial veloc¬ 
ity Vo of each trajectory is sampled from q(v), although 
at later times due to dissipation and driving this distri¬ 
bution will not be conserved. A fraction of this initial 
density contains reactive trajectories. From the survival 
probability of Pr the reaction rate can be expressed as 
the instantaneous flux-over-population. The flux calcu¬ 
lation is formally exact because the DS attached to the 
TS trajectory is recrossing-free. 


A. Harmonic barriers 


When the barrier is harmonic (e = 0), reactive trajec¬ 
tories will cross the moving DS at a time^ 


ft 




■In 


/ Avq - AuAxo \ 

V Avo - AgAxo J ' 


(18) 


The crossing time is a monotonically decreasing function 
of the initial velocity Avg: fast trajectories cross earlier. 
It diverges as Avg AgA^o approaches the stable man¬ 
ifold, and it tends to zero as Avg —?► oo. 

At any time t > 0, the product region Ax > 0, to 
the right of the moving surface, will contain all those 
trajectories that cross the surface at a time < t. These 
are the trajectories that have an initial velocity of at least 


'Cmin = u^(0) -|- Aumin, where t^(Avmin) = t. From this 
condition, we obtain 


AVmin = 


Ke - Age 


Axg. 


(19) 


The population of the product region at time t is there¬ 
fore 


pOC 

Pp(t) = / q{v) 
(^) 


dv, 


( 20 ) 


and the flux across the moving surface is 
dPp 


Fuit) = 


dt 




dv„ 


dt 

dAvm 

dt 


— ArCo (Ag, Ag) 




^gAut — gAsi^2 


( 21 ) 

This result is positive because Axg < 0. 

Alternatively, the flux can be calculated directly from 
the flux integral 

poo 

Fuit) = / dAv Avpt{Ax = 0, Av), ( 22 ) 
Jo 

where pt{Ax,Av) is the density of trajectories in phase 
space at time t. Initially, this density is 

pg{Ax, Av) = S{Ax — Axg) ( 7 (u^( 0 ) -I- Av). (23) 

At later times, it can be obtained from 

Pt{Ax, Av) = pg{Ax[-t), Avfy-t)). (24) 


Here Ax{—t) and Av{—t) denote the phase space point 
reached from Ax, Av by propagating backwards to —t, 
i.e., it is the initial condition that has reached Ax,Av 
at time t. The exponential prefactor accounts for the 
shrinkage of phase space volume: The relative dynamics 
stretches distances at a rate Au in the u direction and by 
a rate Ag < 0 in the s direction. Volumes therefore are 
“stretched” at a constant rate Au + Ag = —7 < 0 , and 
densities must increase accordingly. 

Because the relative dynamics is linear, the equations 
of motion can be solved explicitly. The result is 


Ax{—t) = Qx Ax + Uy Av, 
Av{—t) = bx Ax + by Av, 


with 


Au e - As e 


^-7! - 


Au - As 

AuAg(e-^"*-e-^-*) 
Au ~ As 


67 , — 


g Aui _ g Agt 

A„-A. 

Aue~^"* - Age~^°* 
Au As 


We thus obtain the flux integral 
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pOO 

Fuit) = / dAv Av S{av Av — xp) q{v^{0) + by Av) 

Jo 

= e^* / dAv Av ^ ~ q{vH0) + K Av) 

Jo |ai;| 

= —— — q {v^iO) + by/ayXo) , 

-ay Gy 


(25) 



FIG. 3. The asymptotic product population Pp{oo) of the 
harmonic potential (e = 0) as a function of driving frequency 
n for the fli frequency set. The curves are colored with re¬ 
spect to the value of the initial phase shift: ip = 0 (blue), 
ip = n/2 (cyan), and ip = n (red). For each value of p, the 
dependency of the asymptotic population on the friction pa¬ 
rameter 7 is shown by varying the linestyle: 7 = 0 (solid), 
7 = 1 (dashed), and 7 = 2 (dotted). 


As shown in Fig. O the asymptotic population of the 
product region depends strongly on the frequency of the 
barrier motion fl, the initial phase p, and the friction 7 . 
The asymptotic value is approached according to 

r r’min (^) 

Pp{t) = Pp{oo) — / q{v) dv 

( 00 ) 

= Pp(oo) -b q(Umin(oo)) (Au “ As) Axq 

+0 . (29) 

The rate of approach, i.e., the barrier crossing rate is 

Au - As = y'7^"+4^. (30) 

It depends only on the damping and the shape of the 
barrier, but not on the details of the barrier motion or 
the distribution of initial conditions (unless q{vmin{oo)) 
happens to vanish). 


which can be shown to agree with Eq. m- 

In the limit t ^ 00 , the minimum velocity (1191) is ap¬ 
proximately 

Aumin = As Axo - (Au “ As)Aa:o e“Au-Aa)t 

+0 _ (26) 

As expected, it tends to As Axq, which is the location of 
the stable manifold. Therefore, 

'!^min(oo) = u^(0)-b As Axo = (27) 

The critical velocity is determined by the point of 
intersection between the stable manifold and the line 
a: = xo of initial conditionsThe identification of 
allows the separation of reactive {vq > V^) and non¬ 
reactive trajectories (vo < V^) from initial conditions. 
The stable manifold at t = 0 can be calculated through 
extension of this point to all values of xq and defines a 
critical curve V^. As illustrated in Figs. UKa) and[2Kc), 
V)? is a time-invariant phase space object which separates 
the reactive and nonreactive basins. 

The product population in the long-time limit is 

pCO 

Ppipo) = / q{v)dv. (28) 

Jv'i{0)+\B Axq 


B. Anharmonic barriers 


Analogous to the harmonic case, for an anharmonic 
barrier, we assume that there is a unstable PO with the 
period of the driving. This is similar to the POs used 
by Lehmann et al^^— for the case of thermal activation 
with additive periodic driving. Note that this TS trajec¬ 
tory pl is an exact solution to the equations of motion. 
As pi is an unstable PO, it has stable and unstable man¬ 
ifolds attached. The manifolds are uniquely defined and 
can be calculated perturbatively^l^ or with a numerical 
scheme. The dependence of these manifolds on e and the 
corresponding phase space reaction dynamics are shown 
in Figs. [ 2 Ka) and [ 2 Kc). With increasing anharmonicity, 
also increases due to curvature in the stable mani¬ 
fold. This results in a decrease in fraction of trajectories 
that surmount the barrier leading to products. 

The nonlinear equations of motion (1161) cannot, in gen¬ 
eral, be solved exactly. Let 


^(Fo, to; t) 


/ Pa:{To,to;t)\ 

y7^v(Po; to; t) J 


(31) 


represent the phase space point that is reached at time t 
by a trajectory that starts at Po at time to. Because of 
the external driving, it depends on t and to separately. 
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not only on the difference t — to • The Jacobian matrix of This result is positive because the initial velocity is a 
this trajectory with respect to the initial conditions is decreasing function of the crossing time. 




/ dpx 


dxo 

dvo 

dpv 

dpv 

\dxo 

dvo ) 


(32) 


All derivatives on the right hand side of (15^ are to be 
evaluated at (ro,to;i)- 

Reactive trajectories are those that have an initial ve¬ 
locity Vi > V^, for a critical velocity measured at xq. 
Each reactive trajectory will cross the moving DS at a 
time tci^Vi) and with a velocity Vc- If the crossing time 
decays monotonically from tc{AV^) = oo to tc{oo) = 0 
the inverse function Avi{tc) or Vi{tc) = ul(0) -I- Svi{tc) 
can be obtained. For any crossing time tc > 0, there is 
a unique initial velocity Vi that will lead to a crossing at 
the given time. 

The population of the product region Ax > 0 at time 
tf is therefore, as in Eq. (OUl) , 


pOO 

Ppitc) = / q{v)dv, 


Vi{tc) 

and the flux across the moving surface is 


= -q{vi{tc)) 
= -q{vi{tc)) 


dvj 

dtc 

dAvi 

dtr 


(33) 


(34) 


The flux can also be evaluated directly from the flux 
integral (l 2 ^ 


Fuitc) = / dAv Avpt^{Ax = 0,Av), (35) 


where pt{Ax, Av) is the density of trajectories in phase 
space at time t. Initially, this density is (Eq. (12^ 1 


PoiAx, Av) = S{Ax — Axq) q(v^(0) -I- Av). (36) 


At later times, it can be obtained from Eq. (1241) 


Pt{Ax,Av) = e'^*po{(pa;{Ax,Av,t;0),(py{Ax,Av,t;0)). 

(37) 

Here we have used the general notation for the flow of 
the equation of motion. The exponential accounts for the 
shrinkage of phase space volume and the corresponding 
increase in density. It is the same as in the harmonic 
case: In general, the flow of a differential equation ii = 
f{u) leads to a stretching of volume whose rate is the 
divergence of the vector field f. For Eq. (HU), this rate is 
constant — 7 , so that over time t all volumes will shrink 
by a factor e~^*. 

The flux integral formula (l22l) now reads 


Fuitc) = 



dAv Av J((^a:(0, Av, tc', 0) 


Axo) g(u^(0) -l-v3«(0, Au,tc;0)). 


(38) 


The 8 function requires that the trajectory that reaches 
Ax = 0 , Av at time tc must have started at Axq at time 
0. It produces a single contribution to the integral at 
velocity Avdtc), so that 


Fuitc) = qiv^itti) + Avdtc)) 


AVcitc) 


dipx. 


d^VQ 

tc 


(39) 


where (fviO, Av,tc',0) = Av^tc) and the subscript tc 
indicates that the derivative is to be evaluated at 
(0, Auc(tc), A;0). Similarly, a subscript 0 will be used 
to require evaluation at (Axq, Aviitc),0',tc)- These sub¬ 
scripts indicate derivatives of the flow taken along the 
trajectory from (Axq, Aui(tc)) at t = 0 forward in time 
to ( 0 , Avcitc)) at t = tc (subscript 0 ) and along the same 
trajectory backward in time (subscript tc)- 


To verify that the flux integral (15^ gives the same 
result as (1341) that was obtained from the product popu¬ 
lation, it must be shown that 


dAvi 

Avcitc) 

dtc 

dtpx 



OAvq 

tc 


(40) 


To this end, first note that Av^tc) is defined by the con¬ 
dition 


ifxiAxQ, Aviitc),t)',tc) = 0 . 


Differentiating this condition with respect to tc gives, 


dpx 

dAv 


dAvi 

dtc 



(41) 


0 


0 






















Now difxjdt is the velocity of the trajectory at the end 
point. The second term in Eq. m is therefore Avc{tc)- 
With this result, the condition (l40l) simplihes to 


^ d(px 
dAv f ^ dAv 


(42) 


Under the given assumptions on the geometry, the deriva¬ 
tive on the left hand side is negative: A trajectory that 
arrives at the DS with larger velocity must have started 
further away, i.e., at smaller Ax(0). 

The derivatives occurring in Eq. (|42ll are elements of 
the Jacobian matrices 


J\ta = J(0, AUc(tc),U;0) = 


^ dpx 

dpx 

\ 

dx 

, dv 

tc 

dipv 

dpv 


^ dx 


tc J 


and 


J|o = J{Axo,Avi{tc),0;tc) 


/ difx 

dpx 

\ 

dx 

0 9v 

0 

dify 

dlfiy 



0 

0/ 


respectively. Because these matrices describe variations 
around the same trajectory, taken forwards and back¬ 
wards in time, they are inverse to each other. Eormally, 
this can be shown by taking derivatives of the flow prop¬ 
erty 


#($(ro,0;U),tc;0) = To for all To, 

which says that propagating an arbitrary phase space 
point Fq forward in time by U and back again will return 
the original point. 

By the well known formula for the inverse of a 2 x 2 
matrix, it follows that 



( dify 

difx 

\ 

1 

dv 

0 

0 

det J o 

dify 

dpx 



\ dx 

0 

0/ 


so that 

dipx _ 1 dipx 

dAv det J|o dAv 


C. Dynamics near the TS 


The TS trajectory is a moving saddle point and thus 
trajectories in the neighborhood of Pl can be described 
by a linearization of the equations of motion. In the 
phase space vector relative coordinate AF = (Aa:, Av) 
this linearization is given by 

At = J{t) AT (43) 

where 

( 0 1 

J{t) = 

\ujl + 3e{xi{t) - E{t))‘^ -7 

is the Jacobian of Eq. (HU). The asymptotic decay rate of 
Pn{t) is determined by the behavior of trajectories with 
initial conditions close to the stable manifold. Eor an en¬ 
semble of trajectories constituting an initial phase space 
density po, trajectories that emanate close to (the sta¬ 
ble manifold at t = 0) will persist in the neighborhood 
where (H51) is valid for long times. The decay of these 
trajectories determines the reaction rate. 

The stretching and compression of phase space about a 
PO is known to dictate escape rates in conservative^^— 
and dissipative systems— When J{t) is periodic in sys¬ 
tems of the form of Eq. (l43l) . the rate of deformation 
in the linearized phase space can be quantified through 
calculation of the Floquet exponents— 

To classify the stability of AF^ we consider the dy¬ 
namics of a perturbation vector Acr(t). The equation of 
motion (H51) is linear in Acr{t) and thus it satisfies 

A& = J{t)Acr, Acr(O) =/, (45) 



where I is the 2x2 identity matrix. The principal fun¬ 
damental matrix solution over one period of the driving 
is the monodromy matrix 


^^/AaW(r) Aa(2)(r) 
iACT(i)(r) A(T(2)(r) 


(46) 


A fundamental matrix solution A<j{t) of (I43|) at some 
later time t -I- kT, for k = 1,2,3..., can be obtained as 

A(T{t + kT) = M'^Aait), (47) 


through repeated operation by the monodromy matrix. 

The eigenvalues nis^u of M are the Floquet multipliers. 
The Floquet exponents 


Now 

det J|o = 

is the factor by which phase space volumes shrink during 
time tc- This proves the condition and therefore the 
equality of the two flux formulas. 


^ In (48) 

quantify the stability of AF^ and give the rate of ex¬ 
pansion or contraction of the perturbation of per unit 
time— The TS trajectory has both an unstable /i„ > 0 
and a stable ps < 0 exponent which correspond to 
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stretching and contraction of the initial perturbation in 
the directions of the unstable and stable manifolds, re¬ 
spectively. 

For an arbitrary time interval of length T, trajectories 
that cross the DS in this interval form a strip in the phase 
plane. Trajectories that cross the DS in the next follow¬ 
ing time interval T form a similar strip that that is closer 
to the stable manifold. In the region where the linearized 
system is valid, the phase space density is constant. The 
flux of trajectories through the DS in a given time inter¬ 
val is proportional to the width of the strip that contains 
these trajectories. During sequential periods this width 
decreases by a factor From this it follows 

that, up to periodic modulation, the flux must decay as 
g-lAiu-Ms)* barrier crossing rate is 

= Mm - Ms, (49) 

which expresses the reaction rate in terms of the char¬ 
acteristic Floquet exponents of the TS trajectory. Equa¬ 
tion generalizes Eq. o for the case of an anhar- 
monic barrier. 


V. NUMERICAL RESULTS AND COMPARISON WITH 
THEORY 


The reaction rate of 0 was calculated by simulat¬ 
ing ensembles of n = 10® — 10® trajectories for vari¬ 
ous sets of parameters { 0 , 7 ,e,CT} and following the sur¬ 
vival probability of Pr, as a function of time. A Runge- 
Kutta-Maruyama scheme^S was implemented to perform 
the integration. In the absence of noise {a = 0), this 
algorithm is the well-known fourth-order Runga-Kutta 
method. For all numerical simulations non-dimensional 
parameters were used by choosing units such that the 
barrier frequency Wb and driving amplitude are unity. 
Each trajectory was given an initial position xq = —0.1 
(in the reactant region) and vq was sampled from a Boltz¬ 
mann distribution with ksT = 1. The choice of initial 
conditions is arbitrary as the asymptotic decay rate of 
Pjiit) is independent of the choice of initial distribution, 
suffice that there is enough density about the stable man¬ 
ifold such that a rate exists.— 

The ensemble of n trajectories was evolved through the 
equations of motion ([2]). The normalized reactant popu¬ 
lation was calculated at each time step in the integration 
scheme. An indicator function was employed to follow 
the state evolution of each trajectory. 


hB,[x{t)] 


0, x{t) > x^{t) , 
1, x{t) < x^lt) , 


(50) 


where x^{t) is the configuration space projection of the 
TS trajectory. If for a specific trajectory i, Xi{t) > x^{t) 
that trajectory is in the product state and is not counted 
in the reactant population at time t. The instanta¬ 
neous normalized population of the reactant region can 



t 


FIG. 4. Time dependence of the scaled logarithm of the reac¬ 
tant popnlation, — In [Pnit) — Pr,(oo)], for fli (top), $72 (mid¬ 
dle), and $73 (bottom) with $7/ = 5 for all panels. Values of 
the anharmonic parameter are e G {1, 2, 4, 6, 8,10}. The slope 
of each dashed line is the barrier crossing rate fcf. The color of 
each line corresponds to the respective e value. In all panels, 
parameters are 7=1 and = 0. 


be found by summing over all n trajectories and then 
normalizing by a factor 1 /n, 

1 " 

Pnit) = - ^/iR[a;i(t)]. (51) 

i=l 

Trajectories can only exist in one of two states, reactant 
or product, and so the normalized product population 
Pp = 1 - Pr. 

As shown in Eig. IH the scaled logarithm of the nor¬ 
malized reactant population, — In [PR(t) — Pr(oo)], is ap¬ 
proximately linear in time after an initial transient sec¬ 
tion implying a first-order rate process. The asymptotic 
reaction rate fcf can thus be found as the slope of the 
scaled logarithmic curve in the long-time limit. Periodic 
modulation in the decay of Pnit) was found to become 
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(a) 




(c) 




(e) 


. 4 ^ 




FIG. 5. The barrier crossing rates of systems following the equations of motion m as a function of the anharmonic parameter 
e for various frequency sets fls, driving frequencies fl, and values of friction 7 , as denoted in each panel. The circles denote 
the rates kf calculated from the time evolution of Pnit) through numerically simulation and correspond to the dashed lines 
in Fig. 21 The solid curves are the rates predicted by the difference in the characteristic Floquet exponents fiu — fJ-s of the 
corresponding TS trajectory. 


more prominent for low frequency driving (17/ ^ 2). In 
these cases the global exponential rate was calculated as 
an average over these modulations. 

A comparison between the rates calculated from nu¬ 
merical simulation and rates predicted by Eq. is 
shown in Fig. [S] For all frequency sets fJs and pa¬ 
rameter values, agreement is observed. Underdamped 
(7 < 2 ), overdamped (7 > 2 ), and critically damped 
(7 = 2 ) regimes of a corresponding harmonic well were 
considered. Agreement between the rates persists over 
all ranges of damping. For high frequency driving (17/ > 
uih), the exponential rate can be averaged over several 
periods of driving and modulations in the decay are min¬ 
imal, as illustrated in Fig.|4l Periodic modulations in the 
decay of pR(t) are prominent for low driving frequencies 
(17/ Ri Wb) and the integration of n = 10® trajectories 
resulted in reaching the numerical asymptote Pr,(oo) at 


times less than the period of the external driving. In 
those cases for which the integration time was insufh- 
cient to sample the asymptotic region, a larger number 
of trajectories (n = 10®) were integrated. Each trajec¬ 
tory was integrated to a final time of at least tf = 15 
and, as shown in Fig. 01 Pr(oo) is reached well before 
the end of this sampling window. Increasing the num¬ 
ber of trajectories by an order of magnitude resulted in 
improved convergence of the scaled logarithmic popula¬ 
tion and marginally better agreement between the com¬ 
pared methodologies, as shown in Fig. llKe) for 17 = 1. 
The agreement between the two methods is illustrated 
in Fig.jSjd) for the smaller, non-unity, driving amplitude 
case of 172 . The decreased driving amplitude leads to a 
decrease in the reaction rate. 
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FIG. 6 . The percentage of trajectories that recross the moving dividing surface attached to the DTS trajectory as a function 
of noise strength ct with (a) 7 = 0 . 01 , (b) 7 = 0 . 1 , (c) 7 = 1 , and (d) 7 = 3 for single-frequency (fii) driving and various values 
of e and fl. The black vertical line (solid) marks the noise strength when the fluctuation-dissipation theorem is obeyed. 


VI. CHARACTERIZING NOISY REACTIONS WITH 
THE NOISE-FREE GEOMETRY 

In systems in which the strength of an external driv¬ 
ing force dominates over that of the thermal driving, 
statistical quantities can be approximated by those of 
a corresponding purely deterministically driven system. 
For thermally induced reactions, Lehmann, Reimann, 
Hanggi, and^“— have shown that in the overdamped 
(large- 7 ) regime, when a chemical reaction is forced by 
a periodic field the reaction rate is determined in part 
by the geometry of periodic trajectories in the purely de¬ 
terministic phase space. This work was later extended 
to cases with different scaling behaviors between the 
strength of thermal activation and the strength of the 
external field 

Our goal here is to develop a minimalist theory, appli¬ 
cable at the limit where the magnitude of a noise 
sequence (f) is a small enough perturbation to the pe¬ 
riodic driving i?(t) that the TS trajectory of the noise¬ 
less system (the periodic orbit) gives rise to a DS with 
minimal recrossings. This deterministic TS trajectory 
(DTS trajectory) does not solve the equations of mo¬ 
tion ([5]) with a non-zero value of tr. We therefore distin¬ 
guish the DTS trajectory from the true TS trajectory of 
the noisy system (that we do not compute in this work). 

A principal assumption for the use of the noise-free 
geometry is that the phase space density of the ther¬ 
mal system, and its time-dependence, is approximately 


that of the deterministic system, i.e., pt{Axa, ^Va) ~ 
Pt{Ax, Av). As shown in Fig. [21 for small values of a the 
geometry of the thermal system is similar to that of its 
deterministic counterpart. The rate theory developed in 
Sec. IIV Cl for the deterministic system can therefore be 
applied. This is advantageous in applications such as in 
comparisons with experiments in which the exact noise 
sequence is not known. 

Thermal systems in which the fluctuation-dissipation 
theorem (FDR) is not obeyed due to energy dissipation 
constitute non-equilibrium processes. Formal treatments 
of fluctuation-response in periodically forced systems by 
Teramoto, Harada, and Sas a'^^i'^^ provide insight into the 
rate of energy dissipation in such systems. Green et al^ 
have shown that the rate of energy dissipation is directly 
related to the dynamical entropy of the system. To re¬ 
alize non-equilibrium conditions in the present model re¬ 
action, the damping constant 7 is held constant and the 
strength of the thermal fluctuations cr is increased up 
to the point where the FDR is satished. If the initial 
velocities are drawn from a Boltzmann ensemble with 
k^T = 1 (in dimensionless units), this is the case at 
cr = 7 . If cr < 7 the thermal bath is at a lower tempera¬ 
ture than that of the distribution of initial velocities. In 
this case, the temperature of the ensemble of reactants 
will be continually cooled by the thermal bath. Ther¬ 
mal annealment of the ensemble calls for the develop¬ 
ment of postmodern rate theories which rely singularly 
on geometric properties of phase space, and not the dy- 


















12 
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FIG. 7. Time dependence of the scaled logarithm of the reac¬ 
tant population for systems with single-frequency (fii) peri¬ 
odic and thermal driving for 7 = 0.01 (top), 7 = 0.1 (middle), 
and 7 = 1 (bottom). The color of each line corresponds to a 
specific G value. The decay for systems with various anhar- 
monicites e € {1, 3,10} are shown and denoted in each panel. 
The fundamental driving frequency is 11/ = 5 for all panels. 
For visual clarity, each curve is truncated at a point where 
the data became noisy. 


namics of the ensemble itself. This is in stark contrast 
to the TST assumption of equilibrium distributions in 
metastable states and at the TS. 

The percentage of thermal trajectories that recross the 
DS attached to the DTS trajectory is shown in Fig. [5] for 
varying noise strengths a and constant dissipation rates. 
As shown in Figs. |6}a) and|6}b), a minimal number of 
recrossings occur below and up to the FDR threshold for 
small values of 7. For the 7 = 1 case, shown in Fig. me), 
trajectories persist around the BT for long times, lead¬ 
ing to a larger number of recrossings than observed for 
smaller dissipation rates. For the overdamped dynamics 
(7 = 3), shown in Fig. IHKd), the deterministic DS iden¬ 
tifies reactive trajectories adequately only for weak ther¬ 


mal driving (small cr) and strong anharmonicity. As the 
harmonic limit is approached or in equilibrium systems 
the superimposed DS becomes very poor. 

The decay of the scaled logarithm of the normalized re¬ 
actant population, as calculated with the superimposed 
deterministic DS, is shown in Fig.[7]for various parameter 
values. Over all friction regimes, the population decay of 
the systems with additional thermal driving follows that 
of its deterministic counterpart if the noise strength a 
is sufficiently low. For 7 = 1 , when the strength of the 
thermal driving approaches that of the FDR, a decrease 
in the reaction rate is observed. The data presented in 
Fig. [7] becomes highly oscillatory at long times due to 
recrossings of the DS. For visual clarity each data series 
has been truncated to remove this noisy tail. As observed 
in Figs. |l]and[71 for short times (t ^ 0.3), the decay of 
Pr is non-exponential, signifying temporally global non- 
RRKM kinetic behavior. We obtain the rate from the 
long-time asymptotic decay of of Pr, which is represen¬ 
tative of kinetic experiments in which the concentration 
of a reactant species is directly measured over timei^ 

The thermal rates calculated using the DTS trajectory 
are shown in Fig. [51 As expected by the minimal num¬ 
ber of recrossings shown in Fig. |6l stability analysis of 
the DTS can produce an excellent approximation to the 
rate in thermal environments. Through calculation of 
the error between the numerically calculated rate with 
included noise and the rate given by the Floquet expo¬ 
nents of the DTS trajectory, the extent of applicability 
of the noise-free geometry can be quantified. This error 
is < 3% at 7 = 0.01 over all parameter values. It is < 1% 
for e G {5,10}. Increasing the dissipation by an order of 
magnitude (7 = 0.1) results in the same general trends, 
with all errors generally less than 5%. The exceptions 
occur at the noise strength where the FDR is obeyed 
(cr = 0.1) at e G {1,3} and = 5 for which the error 
« 20%. For 7 = 1.0 and D = 10, all calculated errors are 
less than or on the order of 20%, increasing monotoni- 
cally as a function of cr. As illustrated in Fig. |5}e), at 
lower-frequency driving (II = 5) and large noise (cr = 1), 
the error is between 30% —50%. This suggests a practical 
upper bound to the applicability of the noise-free geom¬ 
etry in estimating the reaction rates in the presence of 
noise. Although not shown, for overdamped dynamics, 
stability analysis of the DTS gives an accurate approx¬ 
imation to the rate only in non-equilibrium small noise 
regimes. 

The calculated errors are on the order of the error 
expected from application of variational transition state 
theory (VTST),^ The presented methodology is advanta¬ 
geous over VTST as it does not require the integration of 
large numbers of trajectories or a flux minimization pro¬ 
cedure. Thus, stability analysis of the DTS trajectory 
offers a simple rate calculation methodology that can be 
readily applied, in weak thermal environments, to driven 
chemical reactions with only prerequisite knowledge of 
the geometry of the energy surface and the functional 
shape of the driving waveform. 
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FIG. 8 . The barrier crossing rates of systems following the equations of motion m as a function noise strength a. The rates 
calculated using the DS attached to the DTS trajectory for single-frequency (fli) driving and various values of e, 7 , and Q 
are shown as circles. The horizontal lines (solid) denote the rates given by the Floquet exponents of the corresponding DTS 
trajectory and are colored according to a respective e value. The black vertical lines (solid) denote the noise strength where 
the fluctuation-dissipation theorem is obeyed. 


VII. CONCLUSIONS 


We have shown that in a model chemical reaction sub¬ 
jected to the influence of forcing from a temporally peri¬ 
odic external field, a recrossing-free dividing surface can 
be constructed over an unstable periodic orbit in the 
region of a moving energetic barrier top. This forms 
the basis for future work on specific driven chemical re¬ 
actions that can be represented by a single collective 
variable for the reaction coordinate under nonequilib¬ 
rium conditions. Potential targets include both substitu¬ 
tion and isomerization reactions in which the governing 
multi-dimensional energy surface can be parameterized 
by a single collective degree of freedom. Other possible 
targets include mechanochemical reactions and stimuli- 
responsive assembly mechanisms when the reaction rate 


is dictated predominantly by geometric properties about 
the moving dividing surface. Generally, force-modified 
and temperature-modulated energy surfaces, deforming 
under the influence of temporally varying forces, moti¬ 
vate the development of rate methodologies that go be¬ 
yond the simplistic equilibrium arguments intrinsic to 
classical equilibrium transition state theory. 

The no-recrossing surface constructed here has been 
shown to persist for strongly anharmonic barriers sub¬ 
jected to single-mode and multi-mode driving waveforms. 
A formally exact rate theory has been developed based on 
the flux of reactive trajectories through this recrossing- 
free surface. It rectifies the principal criterion of tran¬ 
sition state theory for periodically driven chemical reac¬ 
tions. 

To circumvent computationally taxing numerical cal¬ 
culations of the reactive flux through this surface, a rate 
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theory has been developed based on the stability of the 
dividing surface. Strong agreement was observed be¬ 
tween the rate predicted by the Floquet exponents of 
a trajectory defining the phase space evolution on the di¬ 
viding surface, and the rate calculated from simulation of 
a large ensemble of trajectories. Thus, in a periodically 
driven chemical reaction, the asymptotic decay rate of an 
initial distribution of reactants can be extracted directly 
from the stability of the time-varying dividing surface 
irrespective of the dynamics of the reactive population. 

Use of the noise-free geometry to approximate the cor¬ 
responding structure of a driven thermal system has been 
shown to give an excellent approximation to the optimal 
dividing surface if the magnitude of the oscillating force is 
large compared with that from the thermal environment. 
For thermally activated processes, the stability exponents 
of the purely periodically driven system can thus be used 
to predict the reaction rates without an explicit treat¬ 
ment of the thermal dynamics. Extensions of the this 
work to include an explicit treatment of the noise, includ¬ 
ing systems with structured solvents environment a^^i^^ 
and systems displaying fluctuating rates^S are possible 
next steps, and ones which we are currently pursuing. 
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